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Abstract— A variational method is presented for solving eigenvalue problems which arise in connection 
wiih the analysis of convective heat transfer in the thermal entrance region of ducts. Consideration is 
given to both situations where the temperature profile depends upon one cross-sectional co-crdinate 
(e.g. circular tube) or upon two cross-sectional co-ordinates (e.g. rectangular duct). The variational 
method is illustrated and verified by application to laminar heat transfer in a circular tube and a 
parallel-plate channel, and good agreement with existing numerical solutions is attained. Then, appli- 
cation is made to laminar heat transfer in a square duct as a check, an alternate computation for the 
square duct is made using a method indicated by Mi.isaps and Pohlhausen. The variational method 
can, in principle, also be applied to problems in turbulent heat transfer. 


Resume — Une methode variationnelle est presentee pour la resolution des probtemes de valeurs 
propres qui se presentent dans Tanalyse de la transmission de chaleur par convection dans la region 
d'entree des conduites. On considere les deux cas dans lesquels le profil des temperatures depend soit 
d’une seule coordonnce dans la section droite (c’est-a-dire tuyaux circulates) ou de deux coordonnees 
(tuyaux rectangulaires). La methode est verifiee par application a la transmission de chaleur laminaire 
dans un tuyau circulate et dans un canal a faces planes paralteles, on obtient un bon accord avec les 
solutions numeriques existantes. Une application est faite ensuite a la transmission de chaleur laminaire 
dans une conduite carree. A titre de controle, un calcul different pour la conduite carree est effectue 
en utilisant la methode indiquee par Millsaps et Pohlhausen. Cette methode variationnelle peut egale- 
ment etre appliquee aux probldmes de transmission de chaleur turbulente. 


Zusammenfassung— Zur Losung von Eigenwertproblemen des konvektiven Warmeiibergangs beim 
thermischen Einlauf in Kaniile wird eine Variations-methode mitgeteilt. Es wird sowohl der Fall 
betrachtet, dass das Temperatucprofil von einer Querschnittskoordinate (z.B. Kreisrohr), als auch 
jener, bei der es von zwei Querschnittskoordinaten (rechtwinkliger Kanal) abhangt. Die Variations*- 
methode wird erliiutert und auf den laminaren Warmeiibergang im Kreisrohr und im ebenen Spalt 
angewandt, wobei mit den bestehenden numerischen Losungen eine gute Obereinstimmung gefunden 
wird. Sodann wird die Methode auf den laminaren Warmeiibergang im quadratischen Kanal 
angewendet und mit einer Iterationsrechnung nach Millsaps und Pohlhausen verglichen. Prinzipiell 
kann die Variations-methode auch auf Probleme des turbulenten Warmeiibergangs angewendet werden. 


Abstract — JJaeTCH BapnauHOHHbifi mctoa pewemiH 3aAaHii oTbicnaHHH coGctbchh ux 3Ha- 
HeHiitt ypaBHemiH, ocHOBanHufi Ha aHa.iHae KonBeHTHBHoro TemioofiMeHa bo BXOflHoft 
ofijiacTH KanaJioB. PaccMaTpiiBaioTCH ABa cjiyHan: oahh, norAa TeMnepaTypHhitt npo<J)HJib 
3aBHCiiT ot oAHoil HoopAHHaTu nonepeHHoro ceneHHH (HanpiiMep, Kpyrjian Tpyfia), h BTopott 
— ot AByx HoopAHHaT nonepeHHoro ceneHnn (HanpiiMep, nanaji npHwoyrojibHoro npo(J)HJi«). 
BapnauHOHHbift motoa iiJiAiocTpupyeTCH h noATBepwAaeTCH npiiMeHenneM ero h CAynanw 
TenJiooCMeHa npn jiaMHHapHOM noTOKe b npyrAoil Tpyfie h b nJiocKo-napaJiJiejibHOM 
naHaAe; nojiynaeTcn xopouiee coBnaACHne c cymecTByiouuiMH HHCJiennbiMH peuienunMn. 
3aTeM, 3tot MeTOA npiiMeiineTCH k TenJiooOMeHy b ycjioBnnx jiaMimapHoro noTona b naHaae 
c KBaApaTHbiM npotJjM.ieM. B nocAeAHCM CAynae aah kohtpoah npoBeAH pacneT no yna3aHHOMy 
MiuiJicancoM ii Hojibxay3eHOM MeTOAy • B iipiininine BapiiaunoHHbiii mctoa mo>kho TaHHie 
npiiMeniiTb ii k 3aAanaM Ten.iooGMeHa b ycjioBnnx TypGyAenTHoro noTona. 


* At present, Professor of Mechanical Engineering, University of Minnesota Minneapolis 14, Minnesota. 
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NOMENCLATURE 

constants in eigenfunction equation (8); 
half-height of parallel plate channel; 
half-length of side of square duct; 
coefficients in temperature distribution; 
specific heat at constant pressure; 
tube diameter, 2 r 0 ; 

functions of cross-section co-ordinates 
’ in equations (4a) and (10a); 

integrals used in the variational pro- 
’ cedure; 

variational expression, I x — /? 2 / 2 ; 
thermal conductivity; 
normal direction; 
index number of eigenvalue; 

Prandtl number, v/a; 
number of terms in eigenfunction equa- 
tion (8) ; 

heat transfer rate per unit length; 
heat transfer rate per unit area; 
tube Reynolds number, wd/v; 
channel Reynolds number, wa/v; 

/7th eigenfunction; 

functions in eigenfunction equation (8); 
radial co-ordinate; r 0 , tube radius; 
static temperature; T w , wall temperature; 
T b , bulk temperature; T 0 , entering 
temperature ; 

axial velocity distribution; n>, average 
velocity; 

cross-section co-ordinate; X = x/a; 
cross-section co-ordinate; Y = y/a\ 
axial co-ordinate; Z = zfa; 
thermal diffusivity, k/ pc. p ; 
nth eigenvalue; 

, cross-section co-ordinates; 

kinematic viscosity; 
density; 

fully developed temperature distribution 
in square duct; 


Subscript: 

fd , denotes fully developed condition. 


INTRODUCTION 

In a previous paper [1], it was shown that 
variational methods could be successfully applied 
to the computation of the fully developed heat 
transfer characteristics for forced-convection 
flow in passages. Here, attention is directed to 
the thermal entrance region of flow passages. 
Our aim is to formulate and apply a variational 
procedure which may be used to determine 
entrance-region heat-transfer results. To illus- 
trate the method and to establish confidence m 
its predictions, variational calculations are first 
carried out for the circular tube and the parallel- 
plate channel, and comparisons are made with 
exact (numerical) solutions available in the 
literature. Then, the variational method is 
applied to the square duct, for which there are no 
entrance-region calculations in the literature and 
for which an exact solution has not been possible.! 
In all the examples considered, the flow is 
laminar and fully developed, while the heat input 
is uniform along the length of the passage. 
However, the variational method may also be 
applied to the isothermal-wall boundary condi- 
tion and to turbulent flow. 

DESCRIPTION OF THE VARIATIONAL METHOD 

General remarks 

As early as 1885, it was demonstrated by 
Graetz [2] that the temperature solution in the 
thermal entrance region of a passage could be 
formulated as an eigenvalue problem. His 
analysis was concerned only with fully developed 
laminar flow in a circular tube having uniform 
wall temperature. Later, it was found [3-7] that 
the eigenvalue formulation would provide 
entrance-region results for both fully developed 
laminar and turbulent flows in circular tubes and 
parallel-plate channels for either uniform wall 
temperature or uniform wall heat flux. The same 
sort of formulation will apply to other non- 
circular ducts, as will be shown in a succeeding 
section. 


t added in proof: A paper, which appeared after this 
work was completed, gives some approximate results 
for rectangular ducts with boundary conditions different 
from those treated here. This work by S. C. R. Dennis, 
A. McD. Mercer, and G. Poots, appeared in Quart. 
Appl. Math. 17, 285 (1959). 
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To illustrate the way in which eigenvalues 
arise in the entrance-region heat transfer analy- 
sis, consider the problem of laminar flow in a 
circular tube with uniform wall heat flux. As 
shown in [6], the longitudinal variation of the 
wall temperature corresponding to the pre- 
scribed wall heat flux q may be written as 

T w - T 0 __ 11 4z/r o 
qrjk 24 Re Pr 

+ ^ C n K„(r 0 )exp --J 0) 

« = 1 

where T 0 is the temperature of the fluid entering 
the tube and z is the distance measured from the 
entrance of the heated section. In the expression, 
(j* and R„, respectively, represent the eigenvalues 
and eigenfunctions of the following equation: 

r o d I" r dR n "I 

r d (r/r 0 ) |r 0 d (r/r 0 )J 

(2) 

dR n jdr = 0 at r = 0 and r = r 0 . 

Solutions of this homogeneous equation with 
homogeneous boundary conditions are possible 
only for a discrete, though infinite number of /3* 
values. Each fi 2 n which permits a solution is 
called an eigenvalue, and the function R n 
associated with the eigenvalue is called an eigen- 
function. Equation (2) is a special case of the 
general Sturm-Liouville type. From Sturm- 
Liouville theory, it is noted in [6] that: 

JJ [r/r 0 - Wr 0 ) 3 ] [OV'o ) 2 - 

-HrlroY-V24]R n d(rlr 0 ) 
Cn ~" {J [r/r 0 -(r/r 0 )»]/?:rf(r/r5“ 

From this, it is seen that once the eigenfunctions 
R n are known, the coefficients C n of equation (1) 
may be evaluated. So h is apparent that a know- 
ledge of the eigen values and~ eigenfuncti ons _pf 
equation (2) holds the key to the temperature 
solution (I). Since the jB® are all positive, it is 
seen that the series contribution to equation (1) 
dies away for large values of z. Thus, the series is 
only important in the entrance region and 
consequently, the eigenvalue problem is of 


interest only in connection with entrance-region 
heat transfer results. For the uniform wall 
temperature problem, the knowledge of eigen- 
values is not only necessary in the entrance 
region, but also in the fully developed domain.* 

Our goal here is to present and apply a varia- 
tional procedure for solving eigenvalue problems 
which may arise in connection with entrance- 
region heat transfer computations. The varia- 
tional method is especially advantageous for 
the lower eigenvalues, and because of this, it 
serves to compliment the approximate tech- 
niques of [8-10] which work best for higher 
eigenvalues. 

Consideration will be given to two general 
types of eigenvalue problems which may be of 
heat transfer interest. The first type arises whe^- 
the temperature profile depends only on one 
cross-sectional co-ordinate, for example, the 
radial position in a circular tube or the distance 
from the center line of a parallel -pi ate channel. 
We will call this the one-dimensional eigenvalue 
problem. The second type arises when the tem- 
perature profile depends upon two cross-section 
co-ordinates, as in a rectangular duct. We will 
call this the two-dimensional eigenvalue prob- 
lem. Separate presentations and examples will be 
given for each of the two types of problems. 

The one-dimensional eigenvalue problem 

We direct our attention to the Sturm- 
Liouville eigenvalue problem, which will include 
as special cases all heat transfer situations in 
which the temperature depends on only one 
cross-section co-ordinate. The Sturm-Liouville 
problem with which we are concerned is to find 
the eigenvalues! and eigenfunctions R n of the 
following equation: 

U e ^) +fR - +i, ' sK - 0 <4a) 

with the condition that on the boundary surface 
___ R n — 0 or dRJdr) = 0 (4b) 

The functions e, f and Ipnay depend upon the 


* The first eigenvalue, at least, is necessary in the 
fully developed region. 

t Only positive eigenvalues will arise in the problems 
of interest here. 
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independent variable 77. Now, according to the 
calculus of variations (e.g. [11]), there corre- 
sponds to equation (4) the following variational 
expression J : 

J = (5) 

h = S b a [e(dRJdriY -fPQfa (6) 

h = M (7) 

The variational expression J has a very special 
property, namely, that it takes on a stationary 
value, i.e. SJ = 0, when evaluated using a 
and an which satisfy equations (4a) and (4b). 
This characteristic suggests a procedure for 

obtaining approximate solutions for the eigen- 
values and eigenfunctions of equations (4a) and 
(4b) by using the variational expression 7. The 
procedure is based on the Ritz method. Accord- 
ing to this approach, a set of functions R nl , 
R„ 2 , ... is selected, each of which satisfies the 
boundary condition (4b). With these, the nth 
eigenfunction R n is written as 

= A nl R nl + A n2 R, 2 + . . = 2 A ni R vi (8) 

i = 1 

where the A nl , A n2 , ... are constants which 
remain to be determined. As explained later in 
detail, these unknown constants are found by 
imposing the condition that J take on a stationary 
value. An additional condition which must be 
satisfied by the eigenfunctions is that they are 
mutually orthogonal with respect to the weight- 
ing function g, i.e. 

f 0 =j b a gRnRmdr,=0, 

m = n — 1, n — 2, . . 1 (9) 

Thus, for the nth eigenfunction, there are n — 1 
orthogonality equations. As a consequence, it is 
necessary to have at least n terms in the eigen- 
function expression (8) to insure that the varia- 
tional procedure can be fully carried through. 
For each eigenfunction, the variational pro- 
cedure may be reapplied successively using an 
increasing number of terms in the series (8) 
until convergence to a desired accuracy is 
achieved. The convergence is hastened by 
choosing the form of the R ni in accordance with 
any intuition or knowledge one may have about 
the problem. We now turn to describing the 


detailed steps by which the variational procedure 
is applied. 

First eigenvalue: We begin by writing an 
expression for R x in the form of equation (8), 
selecting each function R u in the series to satisfy 
the boundary conditions. With this, the integrals 
A and / 2 , equations (6) and (7), are evaluated, 
and the variational expression J, equation (5), is 
constructed. For the first eigenfunction, the 
orthogonality condition (9) need not be con- 
sidered, since there are as yet no other eigen- 
functions. To find a stationary value of J, the 
expression is differentiated successively with 
respect to each A u and each resulting equation 
is set equal to zero. This provides a set of p 
linear, homogeneous (right-sides equal zero), 
algebraic equations involving the A 1V A 1Z , . . ., 
A lp . A solution to such a homogeneous set is 
possible only if the determinant of the coefficients 
is equated to zero, and from this, there arises 
a polynomial of which fi* is the smallest root. 
Then, returning to the p homogeneous algebraic 
equations, it is possible to find p — 1 of the 
constants A u . The last coefficient is found by 
some other condition such as normalizing, i.e. 
setting / 2 — 1, or else by assigning some arbitrary 
value to R ± at some particular value of rj. As has 
already been noted, the procedure can then be 
repeated using a larger number of terms in the 
series (8) until a desired accuracy is attained. 

Higher eigenvalues: The variational nuthod 
for higher eigenvalues follows a path similar 
to that outlined for the first eigenvalue, except 
that the orthogonality condition (9) must now 
be incorporated into the procedure. This con- 
dition requires that the eigenvalues and eigen- 
functions be found in ascending order; since the 
computation for the /7th eigenfunction R n 
assumes that all preceeding eigenfunctions are 
known. By employing the orthogonality con- 
ditions (9), it is possible to solve for n — 1 of 
the A ni in terms of all the others. Then, turning 
to the expression for /, these particular A ni may 
be eliminated by direct substitution. Thus, the 
number of A ni remaining in J is reduced to 
p ~ (n — 1). The stationary value of J is found 
by differentiating successively with respect to 
each of the remaining A ni and equating each of 
the resulting equations to zero. From these 
p — (n — 1) linear homogeneous equations, the 
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eigenvalue ^ can be found in a manner identical 
to that described in connection with the com- 
putation of the first eigenvalue. Then, the 
constants A ni are computed using these 
p — (n — 1) equations and the orthogonality 
conditions. Because these equations are homo- 
geneous, there will still be one undetermined A ni 
and this can be found by imposing additional 
conditions as previously described. The process 
can then be repeated, taking additional terms 
in the series (8) until a suitable result for jSJ an d 
R n is achieved. 

To help fix these ideas, illustrative examples 
will be given in a later section. 


This characteristic serves as a basis of an 
approximate method for finding the and R n 
which is essentially identical to that which has 
already been described for the one-dimensional 
case. In applying the procedure to the two- 
dimensional case, it is to be remembered that the 
functions R ni which make up the series (8) may 
now depend upon the two co-ordinates ^ and 
773, instead of on the single 77 as before. Also, for 
the two-dimensional situation, the orthogonality 
condition takes the form 

/ 0 = J J gKRm dm dy h = 0, 

area ' .... 

m = n — 1, n — 2, . . .,1 (14) 


The two-dimensional eigenvalue problem 

For the situation where the temperature 
depends upon two cross-sectional co-ordinates, 
the two-dimensional eigenvalue problems which 
lend themselves to solution by the variational 
technique are given by the following equation 
[ 12 ]: 


_0 / a*»\ a / 

d vi \ d vi) + *it\ e dr )2/ 


+//?„ + J8 lgR n = 0 

(10a) 


with the condition that on the bounding surface 

R n = 0 or dRJdN = 0 (10b) 

In this equation, the functions e, f and g may 
depend on both co-ordinates ^ and ij 2 . In a 
manner parallel to the one-dimensional case, 
there exists a variational expression J which 
corresponds to the eigenvalue problem defined 
by equations (10a) and (10b). Rephrasing the 
findings of [12], J may be written as 

J = h~ fik ( 11 ) 


area ( 12 ) 



area 


As before, the variational expression has the 
particular property that it takes on a stationary 
value when evaluated with an eigenvalue and 
eigenfunction of equations (10a) and (10b). 


With these modifications, the detailed directions 
for using the variational method in the one- 
dimensional problem also apply here, and hence 
they need not be repeated. 

An entrance-region heat transfer computation 
for a square duct will be carried out in a later 
section, and this will aid in further clarifying the 
application of the variational procedure to the 
two-dimensional case. 

APPLICATION TO A CIRCULAR TUBE 

As a first illustration of the application of the 
variational procedure, we consider the problem 
of laminar flow in a circular tube with uniform 
wall heat flux. The solution for the wall tem- 
perature as a function of position along the tube 
length has already been given in equation (1), 
while the associated eigenvalue problem is 
defined by equation (2). Now, for the uniform 
heat flux problem, it is well known that the first 
eigenvalue is zero, while R x is a constant 
usually taken as unity. It is also easy to show 
from equation (3) that C x = 0. So, with and 
R x known, we turn our attention to finding the 
second eigenvalue and its corresponding 
eigenfunction /? 2 - 

To apply the variational method, we first 
compare equation (2) with the general form (4) 
and find that 

V - r/r 0 , e =7], f=0, I ~V 2 Y 

With this, the integrals 7 0 , / x and / 2 become 

/o = 1^0 ~v 2 )RnR m dri=0, 

m = n — 1, n — 2, . . 1 (15) 
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h = JJ vidRJdrjy d V , 

7 2 = — i 2 )R 2 n dy ( 16 ) 

The next step is to write R 2 as a sum of terms in 
the form of equation (8). In selecting the func- 
tions R 2i , we take cognizance of the boundary 
conditions (2) and of the fact that the eigen- 
functions can have both positive and negative 
lobes between 77 = 0 and 77 = 1. Functions 
which immediately suggest themselves for the 
role of R 2i are cos 0 7777 = 1, cos 7777, cos 2m i /, etc. 
So as a first approximation for R 2 in the form 
(8), we write : 

R 2 ~ A 21 + A 22 cos 77-77 ( 17 ) 

Then, substituting in equations ( 15 ) and ( 16 ) and 
carrying out the integration yields:* 



(18c) 


Now, using equations (18b) and (18c), we can 
construct the variational expression 


J = h~P\h 

Then, by the orthogonality condition (l 8a), A 21 
can be eliminated from /, giving 



To find the stationary value, we take dJ/dA 22 = 0, 
and from this it follows that 


P\ = 28-997 (20) 

Considering the simplicity of the approximation 
and the relative ease of computation, this result 
is in surprisingly good agreement with the exact 
value [6] : 

PI ~ 25-6796 (21) 

* It is to be remembered that R n-1 = R x = constant. 


The constants A 21 and A 22 remain to be deter- 
mined. From the condition dJ/dA 22 = 0, there 
is obtained an algebraic equation which tells us 
nothing about the constants. But, from the 
orthogonality condition (18a), we have the 
ratio of A 21 to A 22 . With this, the eigenfunction 
expression (17) becomes 

R 2 = j 4 22 (0-087482 + cos 7777) (17a) 

As expected in accordance with what has been 
said in the general presentation of the varia- 
tional method, there remains one undetermined 
constant which must be found from other con- 
ditions. To facilitate comparison with [6], we 
impose the condition that R 2 (0) = 1, and with 
this, equation (17a) becomes 

R 2 = 0-080445 + 0-91956 cos 7777 (17b) 

The value of R at r = r 0 (77 .= 1) plays an 
important role in the wall-temperature com- 
putation, as may be seen from equation (1). 
The first approximation, equation (17b), gives a 
value 7? 2 (1) = —0-83911, as opposed to 
—0-49252 from the exact solution. This compari- 
son strongly suggests that a higher variational 
approximation be carried out for R 2 . 

As a logical refinement of equation (17), we 
add on a term cos 277-77 and write 

^2 — ^21 + ^22 cos 77-77 + A 23 cos 277-77 (22) 

Proceeding as before, 7 0 , 7 X and I 2 are computed 
by integration of equation (22). From the 
orthogonality condition, 7 0 = 0, there is 
obtained 

A 21 = 0-087482 A 22 + 0-303964 A 23 (23) 

Using f x and 7 3 , the variational expression 
J = I 1 — P\I 2 is evaluated to be 

J = 2-46740 ^ 22 - 3-55556 A 22 A 23 + 

' + 9-86960 - >4J 3 - 

- ^ [0-0870045^4 1 2 - 0-0121336 A 22 A 23 + 

+ 0-115501 A \ 3 — 0-25 A\,} (24) 

A 01 can be eliminated from this expression by 
using equation (23), leaving only A 22 and A 23 . 
The stationary value of J is achieved by taking 
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8J/8A 2 2 = 0, dJ/8A 23 = 0. From this, there 
results: 


= 0 (4-93480 - 0- 1 70 1 82 p l)A 22 + 

(/yi 22 

+ (-3-55556 + 0-0254293 13 *)A 23 = 0 


dj 

3^23 


= 0 -* (-3-55556 + 


r (25) 


+ 0-0254293/8 \)A 22 + (19-7392 - 
— 0 184805 ^)^ 23 = 0- 


A nontrivial solution of this pair of linear, 
homogeneous equations is possible only if 


(4-93480 - 0-170182 $) 

(-3-55556 + 0-0254293 p\) _ 

(—3-55556 + 0-0254293 /3jj) 

(19-7392 - 0-184805 ]SJ) 

This determinant yields a quadratic equation for 
PI, the smallest root* of which is 


PI = 25-6956 (26) 


This is in remarkably close agreement with the 
value 25-6796 from the exact solution. Now, 
turning to the determination of the constants 
A 2i , there is first obtained from either of equa- 
tions (25) the ratio A 2 JA 22 . Next, from equation 
(23), the ratio A 21 /A 22 is found. Finally there is 
imposed the condition that R 2 ( 0) = 1. With this 
information, all three constants may be calcu- 
lated and the expression for R 2 becomes 


r 2 = 0-109207 + 0-746309 cos 77-17 + 

+ 0-144484 cos 2777 ? (22a) 


Utilizing equation (22a) for R 2 , the coefficient. C 2 
is computed as 0*40680, which deviates by only 
0*8 per cent from the value 0-40348 given by the 
exact solution. 

This illustration demonstrates that the varia- 
tional procedure is capable of providing excellent 
heat transfer results. The high level of agreement 
indicates that there is no need to refine R 2 
further. However, if an exact solution had not 
been available for comparison purposes, it 
would have been necessary to take a higher 
approximation for R 2 to establish the level of 
accuracy. 

If desired, the variational method could now 
be applied to the computation of the next eigen- 
function i? 3 . No essential changes in the method 
are necessary, except that the expression for 
i ? 3 in the form of equation ( 8 ) would likely 
contain additional cosine terms. Also, in the 
orthogonality condition (15), we would use 
equation (22a) for the known eigenfunction R 2 , 
while R 1 is a constant. However, since our 
purpose in considering the circular tube has 
only been to illustrate the variational' method, 
the computation for the higher eigenvalues will 
not be carried out. 

APPLICATION TO A PARALLEL PLATE CHANNEL 

As a second example which may serve to 
establish confidence in the variational procedure, 
we turn to the problem of laminar flow in a 
parallel-plate channel with uniform wall-heat 
flux. Corresponding to the prescribed heat 
flux, the variation of the wall temperature along 
the length is given by: 


The value of R 2 at the tube wall (77 = 1) has 
been already noted as playing an important role 
in the heat transfer results. Equation (22a) 
gives a value 7? 2 (1) = —0*49262, which is in 
excellent agreement with the result —0-49252 
from the exact solution. In addition to R n (\), the 
solution for the wall temperature variation as 
-given by equation ( 1 ) depends^ on the coefficients 
C n . These constants can be computed' from the 
quotient of integrals given in equation (3). 


* Since the value of associated with the exact 
solution is an absolute minimum, it is clear that the best 
result is achieved by selecting the smallest root. 


Tw ~ To _ I? , zja ' 
qajk 35 Re a Pr 

CO 

+ zL c " RM exp { ~ (3/2 Yk^r ^ (27) 
«= 1 

The eigenvalues R n and eigenfunctions p n are 
"found from* the following-homogeneous system _ 

++ 2 + RnPl [1 - (y/af] = 0 (28a) 

d(y/a) 2 

dR„/dy = 0 at y = 0 and y — a (28b) 
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while the coefficients C n are computed from the 
ratio: 

JS [1 -(yl°) 2 ] [fO'/tf) 2 - Ky/af - 39/280 ]R n d(y/a) 
ft [1 -(yla?]Rld(yla) 

(29) 

Here again is a problem ideally suited for the 
variational procedure. Comparing equation 
(28a) with the general one-dimensional form 
(4a), it is seen that 

rj=yla, e = 1, /= 0, g = \ — t] 2 (30) 

and the integrals / 0 , j\ and / 2 become 

m = n — 1, n — 2, . . 1 (31a) 

h = JJ ( dRJdr)) 2 dr), h = Jo (1 — ^? 2 )^ 2 dr) (31b) 

Since we are dealing with the case of uniform 
wall heat flux, it follows that = 0, R ± — con- 
stant, and Q = 0. So, attention can immedi- 
ately be directed to the second eigenvalue R 2 . 
Drawing, upon our experience with the circular 
tube, we start out by writing R 2 in the form of 
equation (22). Following through the operations 
as before, it is found that 

PI = 18-39, i? 2 = -0*3157 + 

+ 1*124 cos 7T7) + 0-1924 cos Ittk) (32) 

The numerical solution of [7] yields a value of 
18-38 for p\. For 7? 2 (1), which is needed in the 
wall-temperature computation, [7] gives —1-27 
as compared to —1-25 from the variational 
solution (32). No results for C 2 are provided in 
[7] and so this computation is omitted here. 

From the excellent level of agreement which 
has been demonstrated, one can draw a real 
feeling of confidence in the utility of the varia- 
tional procedure. 

APPLICATION TO A SQUARE DUCT 
For the previous examples which have been 
presented, there exist numerical solutions in the 
literature. These illustrations have been useful 
in establishing a feel for the accuracy which 
could be achieved by the variational method. 
Now, we turn to a situation for which there are 
no entrance-region calculations in the literature. 
Consideration is given to a laminar flow in a 


square duct, Fig. 1, having a heat flux which is 
uniform both along the length and around the 
periphery. Since this physical problem is noL 
treated elsewhere, we start from first principles. 



Fig. 1 . Co-ordinate system for square duct. 

The energy equation appropriate to the fully 
developed laminar flow of an incompressible, 
constant property fluid is 


P c v w 


dT 

dz 


/ d 2 T 3 2 T\ 

\w + ~df) 


(33) 


where axial heat conduction has been neglected 
compared to transverse conduction. Far down 
the duct, dT/dz becomes a constant and the 
condition of fully developed heat transfer is 
achieved. A solution for the temperature distri- 
bution in the fully developed regime has been 
found by variational means in [1] as follows. 


T-T 0 2 Z 

Q/Sk ~ Re a Pr + ^ ^ 34a ) 

x 2 4 - Y 2 

</> = j -0-16032 [(X 2 - l) 2 + 

+ (Yl- l) 2 ] -0-13290 

[(X 2 - 1) (Y 2 - l)] 2 + 0-06803 (34b) 


where T 0 is the temperature of the fluid entering 
the duct, Q is the heat transfer per unit length 
and X, Y and Z are dimensionless co-ordinates. 

Now, turning to the entrance-region, we 
propose a solution in the following form which 
will apply anywhere throughout the duct: 


T-T 0 _ 2Z , ^ , 

2/8/c Re a Pr 9 ^ 

00 

+ 2 C ’ S ’ eXP (-^Fr Z ) < 35 > 

>»= 1 
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In this equation, R n depends on both cross-sec- 
tion co-ordinates X and Y. To find the governing 
equation for R n , (35) is introduced into the 
energy equation (33), from which it folJows that: 

U - + = 0 (36a) 

It is also easy to show that: 

dR/dX = 0 at X± 1, 

dR/dY = Oat Y=± 1 (36b) 

So, once again we have an eigenvalue problem. 
The first thought about solving equation (36) is 
to try separation of variables, i.e. to suppose that 
R n is a product of a function of X and a function 
of Y. However, this approach is not successful 
since the velocity distribution w/w is not in a 
product form. Hence, it is necessary to deal with 
equation (36a) as it stands. 

The eigenvalue problem as represented by 
equations (36a) and (36b) is well-suited to be 
attacked by the variational procedure. Com- 
paring equation (36a) with the general two- 
dimensional form (10a), it is seen that: 

Vi = V 2 = Y, e = 1, /= 0, g = w/w (37) 

With this, the integrals / 0 , and I 2 as given by 
equations (14), (12) and (13) become 

/o=4 JJJS (w/w)R n R m dXdY = 0, 

m — n — 1 , n — 2, . . . , 1 (38a) 

h = 4 JJ JJ [(dR n /dXy+(dR n ldYy]dXdY (38b) 

= 4 ^{w/w)RldXdY (38c) 

Since we are dealing with the uniform heat flux 
case, the first eigenvalue PI is zero, while 
R x = constant and C x = 0. So consideration 
can be immediately given to the second eigen- 
value R 2 . In selecting a trial function for R 2 , we 
are guided by prior experience with the parallel- 
plate channel. For that configuration, an eigen- 
value expression in the form of equation (22) 
proved to “be^quite satisfactory. Generalizing- 
equation (22) to the two-dimensional case, we 
write 

R 2 = 4 2l + A 22 cos t tX cos ttY + 

+ ^23 [cos I** cos 2 t tY + cos 2ttX cos t t Y] (39) 


where the constants A 21 , A 22 and A 23 remain to be 
determined by use of the variational method. 
The first step is to compute the I integrals of 
equations (38). To carry out the integrations, a 
knowledge of the velocity distribution w/w is 
required. An accurate velocity profile has been 
determined by variational means [1] as follows: 

w/w = (X 2 - 1) (Y 2 - 1) [2-0983 + 

+ 0-29181 (X 2 + Y 2 ) + 0-87546 X 2 Y 2 ] (40) 

Utilizing this expression in conjunction with 
equation (39), the / integrals are found to be 

/ 0 = A 21 + 0-082665 A 22 - 

— 0-0432577 A 23 ~ 0 (41a) 

h = 2^ [Ah + 5A\ Z ] (41b) 

I 2 = -4A 2l +’ 0-83989 A\ % + 

+ 2-0001 A\ 3 + 1T6410 A 22 A 23 (41c) 

Then, the variational expression J = — f3 2 f 2 

is constructed. Utilizing equation (41a), it is 
possible to eliminate A 21 from /, leaving only 
A 22 and A 23 . The stationary value is found by 
setting 3J/dA 22 = 0 and dJ/dA 23 = 0. From this, 
there is obtained 

= 0 -»■ (39-478 - 1-6251 ft )A a - 

O A 22 

— (1-1927 Pl)A 23 — 0 (42a) 
~ =0-» (-1-1927^2 + 

0^23 

+ (197-39 - 3-9852 p\)A 23 = 0 (42b) 

By setting the determinant of the coefficients 
equal to zero, there results a polynomial in 
the smallest root of which 

PI = 20*929 (43) 

provides the second eigenvalue. Next, the 
constants A 21 , A 22 -and A 23 are ^ found by utilizing 
either of equations (42) in conjunction with 
(41a) plus an additional condition. In this in- 
stance, as a matter of variety, we will impose the 
condition that the eigenfunction be normalized, 
i.e. I 2 = 1. There is thus sufficient information 
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to determine all the constants, and the final ex- 
pression for the eigenfunction R 2 becomes: 

R 2 = 0 067685 - 0-92477 cos nX cos 77 Y- 
- 0-20252 [cos 77^ cos 2 t tY + 

+ cos 2ttX cos 77 Y] (44) 

The last bit of information needed to compute 
temperature distribution results from equation 
(35) is the constant C 2 . Imposing the condition 
that T = T 0 at the entrance section (Z = 0), 
equation (35) becomes 

l C n R n = - * (45a) 

n= 1 


is prescribed, the information of greatest 
interest is the resulting wall temperature. In this 
instance, where the heat flux is everywhere 
uniform, the wall temperature will vary around 
the periphery in any cross-section, as well as 
along the length. Because of the symmetry of 
the problem, consideration need only be given to a 
typical part of the wall: X = 1, 0 < Y ^ 1. 
The wall-temperature result can be put in a 
convenient form by first introducing the bulk 
temperature T b : 


Qj%k Re a Pr 


Then, multiplying through by ( w/w)R m and 
integrating over X and Y from 0 to 1, it follows 
by use of the orthogonality condition (38a) that 


JoJo i w lw)R n <j> dX dY 
J&ftO *tmi dX dY 


(45b) 


and then, it may be observed that under fully 
developed conditions 


(T-T b ) fd _ A 


(48) 


With these, equation (35) can be rephrased as: 


Carrying out the integration for R 2 , it is found 
that 

C 2 = —0-2467 (46) 

This completes the computation for the second 
eigenfunction. To check the level of accuracy of 
the results, the variational procedure might be 
reapplied to an expression for R 2 which contains 
additional terms. But, we have decided instead 
to redo the problem using a completely different 
approximation procedure. The method is a 
modification of an idea presented by Millsaps 
and Pohlhausen [13] and its application to the 
current problem is discussed briefly in the 
Appendix. The results of this alternate compu- 
tation are in very good agreement with those of 
the variational method. 

With the second eigenvalue at our disposal, 
we can now turn to a discussion of the heat 
transfer results. The temperature distribution 
corresponding to the prescribed heat flux is 
given by equation (35), where </> represents the 
fully developed solution as written in equation 
(34), while C 2 , R 2 and are given by equations 
(46), (44) and (43), respectively. Since we have 
only one term of the series (C x = 0), attention 
must be directed to that portion of the duct near 
the fully developed region. When the heat flux 


T-T b 
(T — T b ) fd 


= 1 + 


2 CnRn CXP ( RelPr 2 ) 

n — 1 

* 


(49) 


Finally, at the wall (X = 1, 0 < Y < 1), there 
is obtained 



where the series has been truncated after n = 2. 
The group C 2 {R 2 j<t >) A = 1 as evaluated from the 
variational method has been plotted as a solid 
line on Fig. 2. In the region near the mid-wall 
(Y = 0), the group has negative values; while 
near the corner, the group has positive values. 
Using this informat on in conjunction with 
equation (50), it is seen that for locations near 
the comer, the wall-to-bulk temperature differ- 
ence is greater in some part of the entrance- 
region than it is in the fully developed region. 
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Of course, without additional terms in the 
series, it is not possible to know the extent of the 
entrance region where ( T w - T b ) > (T w - T b ) fd . 
This finding is somewhat surprising in view 
of previous experience with one-dimensional 
geometries such as the circular-tube and parallel- 
plate channel. In those configurations, where 
peripheral variations are absent, T w — T b in the 
entrance-region is always less than the fully 
developed temperature difference. In the present 
instance, where there are peripheral variations, 
it is the feeling of the writers that T w — T b is 
not the governing driving force for heat transfer 
at local positions around the periphery. Hence, 
there would be no reason to expect that 
the longitudinal variation of T w — T b should 
follow the same trend here as was found in the 
simpler geometries. The fact that the longitudinal 
temperature variations are different in different 
parts of the cross-section suggests that there is 
little utility in computing average Nusselt 
numbers. 

The findings "derived ~ fronr the variational - 
method are closely confirmed by the alternate 
computation of the Appendix, as may be seen 
from the dashed line of Fig. 2. This good agree- 
ment increases our confidence in the results, but 
ultimate confirmation awaits some more exact 


solution, an approach to which is currently not 
apparent to the writers. 

CONCLUDING REMARKS 
The examples considered here are meant to 
illustrate the method and by no means exhaust 
the problems to which the variational technique 
can be applied. 

Although the illustrations were concerned 
with the uniform heat flux case, solutions for 
uniform wall temperature are obtained with 
equal facility. In the latter case, the functions 
R ni which make up the eigenfunction expression 
(8) would be selected to have zero values at the 
wall, rather than zero derivatives as in the former. 
It is also worth noting that the first eigenvalue 
for the uniform wall temperature situation will 
not be zero, nor will R x be a constant. Aside 
from these matters, the method is applied exactly 
as has been illustrated. 

In principle, the variational procedure could 
be applied to problems in turbulent heat transfer. 
However^ no^ work has thus far been done along 
these lines. 

APPENDIX 

Alternate Eigenvalue Calculation 
According to the idea of Millsaps and Pohl- 
hausen [13], the eigenfunctions for a laminar, 
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convective, heat transfer problem are expanded 
in terms of the eigenfunctions of the slug-flow 
^£l>lem. For the square duct, such an expan- 
sion for R 2 which is truncated after three terms is 



R 2 = D x cos ttX cos t tY + 



D 2 ( COS ttX COS 277 Y + COS 2nX COS 77 Y) + 

+ Z) 3 ( cos 2 t 7* cos 2t t Y) (51) 


where D l9 D 2 and Z) 3 remain to be determined. 
Introducing this expression into the governing 
equation (36a) gives : 


2tt 2 D x COS ttX COS 77 Y + 

+ 57t 2 D 2 (cos nX cos 2tt Y -j~ 

+ cos 277* cos 77 Y) + 
+8t7 2 D 3 (cos2t7*cos 2t7 Y) = £f(w/vv)/? 2 


Then, equation (52) is multiplied through by 
cos 77* cos 77 y and integrated over * and * 
from 0 to 1. Next, equation (52) is multiplied by 
(cos 77* cos 277*+ cos 277* cos v *) and inte- 
grated over the same range. Finally, this same 
procedure is carried through with cos 277* 
cos 2 t 7 Y. The result of these operations is three 
linear, homogeneous, algebraic equations for 
the D l9 D 2 and Z) 3 , the coefficients of which 
contain £ 2 . To obtain a nontrivial solution, the 
determinant of the coefficients is equated to 
zero, and this gives: 


j8* = 20-222 (53) 

which is rather good agreement with equation 
(43) from the variational procedure. The 
D-values are then found by returning to the 
homogeneous algebraic equations, from which 
two of the three D ' s can be determined in terms 
of the third one. The third D value may be left 
unspecified or else arbitrarily assigned. This will 
have no effect on the final result for the desired 
quantity C n R„ (see equation (35)), since a 
change in the level of R n will be automatically 
compensated by a corresponding change in C n . 
Then, turning to the computation of C 2 , the R 2 
expression (with D-values now known) is intro- 
duced into equation (45b) and the integration 
carried out. With this, the final result for the 
C 2 R 2 product is 


C 2 R 2 = 0-2140 cos 77* cos 77 * + 

+0 01 740 cos 277* cos 277 * + 

+0-04698(cos 77* cos 2nY + 

+ cos 2 t7* cos 77 Y) (54) 

Since the form of (54) is somewhat different 
from that of the variational eigenfunction, 
comparisons are best made by evaluating the 
expressions at specific *- and *-values. Fig. 2 
shows a comparison made along the wall, 
* = 1, 0 < Y < 1, with quite good agreement. 
A similar comparison has been made along the 
center-line in the fluid, * = 0, 0 < * < 1, with 
about the same level of agreement. 
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